Template author: Morgan Carr-Markell
Template last modified on: Feb 28, 2023
Notebook modified by: Katie Wang
Notebook last modified on: Mar 1, 2023
The code and many descriptions in this notebook come from a tutorial posted on the (National Ecological Observatory Network) NEON website and originally from the Data Carpentries. You can find the original tutorial here: https://www.neonscience.org/resources/learning-hub/tutorials/introduction-working-raster-data-r#toggle-26
Important note: Remember to:
Otherwise, you will have no data to work with!
Before you start the tutorial on loading, plotting, and processing raster data in R, go to this NEON tutorial and read just to the end of the section called “Gridded, or Raster, LiDAR Data Products”
What do DTM, DSM, and CHM stand for?
DTM: Digital Terrain Model.
DSM: Digital Surface Model.
CHM = DSM - DTM; it’s the Canopy Height Model.
What are the two main formats of LiDAR data files and how are
they different?
One is the .las file format for LIDAR point clouds. These are
collections of points with (x, y, z) values and other attributes, and
are not aligned with a particular grid. The other is the gridded or
raster data format, which is a grid of cells all of the same size. Each
cell represents a certain amount of area on the ground.
Original Authors: Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul
Original Last Updated: Apr 8, 2021
In this tutorial, we will review the fundamental principles, packages and metadata/raster attributes that are needed to work with raster data in R. We discuss the three core metadata elements that we need to understand to work with rasters in R: CRS, extent and resolution. It also explores missing and bad data values as stored in a raster and how R handles these elements. Finally, it introduces the GeoTiff file format.
options("rgdal_show_exportToProj4_warnings"="none") # to suppress warning messages
library(raster)
Loading required package: sp
library(rgdal)
Please note that rgdal will be retired during 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-4, (SVN revision 1196)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
Path to GDAL shared files: C:/Users/achro/AppData/Local/R/win-library/4.2/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
Path to PROJ shared files: C:/Users/achro/AppData/Local/R/win-library/4.2/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-0
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(rasterVis)
Loading required package: lattice
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0 ✔ purrr 1.0.1
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.5.0
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::select() masks raster::select()
# Load raster into R
DSM_HARV <- raster("NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
# View raster structure
DSM_HARV
class : RasterLayer
dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
source : HARV_dsmCrop.tif
names : HARV_dsmCrop
values : 305.07, 416.07 (min, max)
# plot raster
# note \n in the title forces a line break in the title
plot(DSM_HARV,
main="NEON Digital Surface Model\nHarvard Forest")
Raster data can be continuous or categorical. Continuous rasters can have a range of quantitative values. Some examples of continuous rasters include:
The raster we loaded and plotted earlier was a digital surface model, or a map of the elevation for Harvard Forest derived from the NEON AOP LiDAR sensor. Elevation is represented as a continuous numeric variable in this map. The legend shows the continuous range of values in the data from around 300 to 420 meters.
Some rasters contain categorical data where each pixel represents a discrete class such as a landcover type (e.g., “forest” or “grassland”) rather than a continuous value such as elevation or temperature. Some examples of classified maps include:
Note: Skipped section on creating categorical plot because it is covered in more detail below
Raster data can come in many different formats. In this tutorial, we will use the geotiff format which has the extension .tif. A .tif file stores metadata or attributes about the file as embedded tif tags. For instance, your camera might store a tag that describes the make and model of the camera or the date the photo was taken when it saves a .tif. A GeoTIFF is a standard .tif image format with additional spatial (georeferencing) information embedded in the file as tags. These tags can include the following raster metadata:
In this tutorial we will discuss all of these metadata tags.
More about the .tif format:
Note: I am skipping the CRS part of the tutorial because we covered that last week, but if you want a refresher you can look at that section of the online tutorial here
The spatial extent is the geographic area that the raster data covers.
The spatial extent of an R spatial object represents the geographic “edge” or location that is the furthest north, south, east and west. In other words, extent represents the overall geographic coverage of the spatial object.
A raster has horizontal (x and y) resolution. This resolution represents the area on the ground that each pixel covers. The units for our data are in meters. Given our data resolution is 1 x 1, this means that each pixel represents a 1 x 1 meter area on the ground.
The best way to view resolution units is to look at the coordinate reference system string crs(). Notice our data contains: +units=m.
Note: You’ll see a warning about the Proj.4 representation being deprecated. You can ignore that.
crs(DSM_HARV)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16018]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Note: I am skipping the raster min and max values section
Raster data often has a NoDataValue associated with it. This is a value assigned to pixels where data are missing or no data were collected.
By default the shape of a raster is always square or rectangular. So if we have a dataset that has a shape that isn’t square or rectangular, some pixels at the edge of the raster will have NoDataValues. This often happens when the data were collected by an airplane which only flew over some part of a defined region.
In the image below, the pixels that are black have NoDataValues. The camera did not collect data in these areas.
# stack() is a raster function to create a rasterStack object
RGB_stack <- stack("NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
# Create an RGB image from the raster stack
par(col.axis = "white", col.lab = "white", tck = 0)
# Note: col.axis and col.lab set the colors for the axis annotation and x/y labels
# tck is the length of tick marks on the axes
# The code below plots an RGB color image using layers in a rasterStack object
# We have to tell the function which layers in the stack contain red green and blue
plotRGB(RGB_stack,
r = 1,
g = 2,
b = 3,
axes = TRUE,
main = "Raster With NoData Values\nRendered in Black")
Note: Skipped section on setting NoData values to NA
The assigned NoDataValue varies across disciplines; -9999 is a common value used in both the remote sensing field and the atmospheric fields. It is also the standard used by the National Ecological Observatory Network (NEON).
If we are lucky, our GeoTIFF file has a tag that tells us what is the NoDataValue. If we are less lucky, we can find that information in the raster’s metadata. If a NoDataValue was stored in the GeoTIFF tag, when R opens up the raster, it will assign each instance of the value to NA. Values of NA will be ignored by R as demonstrated above.
Bad data values are different from NoDataValues. Bad data values are values that fall outside of the applicable range of a dataset.
Examples of Bad Data Values:
Sometimes a raster’s metadata will tell us the range of expected values for a raster. Values outside of this range are suspect and we need to consider than when we analyze the data. Sometimes, we need to use some common sense and scientific insight as we examine the data - just as we would for field data to identify questionable values.
We can explore the distribution of values contained within our raster using the hist() function which produces a histogram. Histograms are often useful in identifying outliers and bad data values in our raster data.
# view histogram of data
hist(DSM_HARV,
main = "Distribution of Digital Surface Model Values\n Histogram Default: 100,000 pixels\n NEON Harvard Forest",
xlab = "DSM Elevation Value (m)",
ylab = "Frequency",
col = "wheat")
Warning: 4% of the raster cells were used. 100000 values used.
Notice that a warning is shown when R creates the histogram.
This warning is caused by the default maximum pixels value of 100,000 associated with the hist function. This maximum value is to ensure processing efficiency as our data become larger!
We can define the max pixels to ensure that all pixel values are included in the histogram. USE THIS WITH CAUTION as forcing R to plot all pixel values in a histogram can be problematic when dealing with very large datasets.
# View the total number of pixels (cells) in is our raster
totalCells <- ncell(DSM_HARV)
Now we can add that as a parameter, maxpixels:
# create histogram that includes with all pixel values in the raster
hist(DSM_HARV,
maxpixels = totalCells,
main = "Distribution of DSM Values\n All Pixel Values Included\n NEON Harvard Forest Field Site",
xlab = "DSM Elevation Value (m)",
ylab = "Frequency",
col = "wheat4")
Note that the shape of both histograms looks similar to the previous one that was created using a representative 10,000 pixel subset of our raster data. The distribution of elevation values for our Digital Surface Model (DSM) looks reasonable. It is likely there are no bad data values in this particular raster.
The Digital Surface Model object (DSM_HARV) that we’ve been working with is a single band raster. This means that there is only one dataset stored in the raster: surface elevation in meters for one time period.
A raster dataset can contain one or more bands. We can use the raster() function to import one single band from a single OR multi-band raster. We can view the number of bands in a raster using the nlayers() function.
# view number of bands
nlayers(DSM_HARV)
[1] 1
However, raster data can also be multi-band meaning that one raster file contains data for more than one variable or time period for each cell. By default the raster() function only imports the first band in a raster regardless of whether it has one or more bands.
Remember that a GeoTIFF contains a set of embedded tags that contain metadata about the raster. So far, we’ve explored raster metadata after importing it in R. However, we can use the GDALinfo(“path-to-raster-here”) function to view raster metadata before we open a file in R.
Note: In this notebook, you’ll need to click through the three data frames to see all of the info about this rasterLayer.
# view attributes before opening file
GDALinfo("NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif")
Warning: GDAL support is provided by the sf and terra packages among othersWarning: GDAL support is provided by the sf and terra packages among othersWarning: GDAL support is provided by the sf and terra packages among others
rows 1367
columns 1697
bands 1
lower left origin.x 731453
lower left origin.y 4712471
res.x 1
res.y 1
ysign -1
oblique.x 0
oblique.y 0
driver GTiff
projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
file NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_dsmCrop.tif
apparent band summary:
apparent band statistics:
Metadata:
AREA_OR_POINT=Area
Notice a few things in the output:
It is ideal to use GDALinfo to explore your file before reading it into R.
Note: rgdal is being replaced by functions in sf and terra/stars so this won’t be true in the future.
GDALinfo("NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
Warning: GDAL support is provided by the sf and terra packages among othersWarning: GDAL support is provided by the sf and terra packages among othersWarning: GDAL support is provided by the sf and terra packages among others
rows 1367
columns 1697
bands 1
lower left origin.x 731453
lower left origin.y 4712471
res.x 1
res.y 1
ysign -1
oblique.x 0
oblique.y 0
driver GTiff
projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
file NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif
apparent band summary:
apparent band statistics:
Metadata:
AREA_OR_POINT=Area
Does this file have the same CRS as DSM_HARV?
Yes; they both have proj=utm.
What is the NoDataValue?
NoDataValue is -9999.
What is resolution of the raster data?
The x and y resolution are both 1.
How large would a 1x1 pixel area would be on the Earth’s
surface?
Because res.x and res.y are 1, and the units of the projection system
are meters, each pixel is a 1m x 1m area on Earth.
Is the file a multi- or single-band raster?
This file has 1 band.
Notice: this file is a hillshade, which uses information about elevation to add shadows to make changes in elevation clearer to a map viewer.
Original Authors: Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul
Original Last Updated: Apr 8, 2021
This tutorial reviews how to plot a raster in R using the plot() function. It also covers how to layer a raster on top of a hillshade to produce an eloquent map.
In this tutorial, we will plot the Digital Surface Model (DSM) raster for the NEON Harvard Forest Field Site. We will use the hist() function as a tool to explore raster values. And render categorical plots, using the breaks argument to get bins that are meaningful representations of our data.
First, let’s plot our Digital Surface Model object (DSM_HARV) using the plot() function. We add a title using the argument main=“title”.
# Plot raster object
plot(DSM_HARV,
main="Digital Surface Model\nNEON Harvard Forest Field Site")
We can view our data “symbolized” or colored according to ranges of values rather than using a continuous color ramp. This is comparable to a “classified” map. However, to assign breaks, it is useful to first explore the distribution of the data using a histogram. The breaks argument in the hist() function tells R to use fewer or more breaks or bins.
If we name the histogram, we can also view counts for each bin and assigned break values.
# Plot distribution of raster values
DSMhist <- hist(DSM_HARV,
breaks = 3,
main = "Histogram Digital Surface Model\n NEON Harvard Forest Field Site",
col = "wheat3", # changes bin color
xlab = "Elevation (m)") # label the x-axis
Warning: 4% of the raster cells were used. 100000 values used.
Warning message!? Remember, the default for the histogram is to include only a subset of 100,000 values. We could force it to show all the pixel values or we can use the histogram as is and figure that the sample of 100,000 values represents our data well.
# Where are the breaks and how many pixels in each category?
DSMhist$breaks
[1] 300 350 400 450
DSMhist$counts
[1] 31676 67862 462
Looking at our histogram, R has binned out the data as follows:
300-350m, 350-400m, 400-450m
We could specify different breaks, if we wished to have a different distribution of pixels in each bin.
We can use those bins to plot our raster data. We will use the terrain.colors() function to create a palette of 3 colors to use in our plot.
The breaks argument allows us to add breaks. To specify where the breaks occur, we use the following syntax: breaks = c(value1, value2, value3). We can include as few or many breaks as we’d like.
# plot using breaks.
plot(DSM_HARV,
breaks = c(300, 350, 400, 450),
col = terrain.colors(3),
main = "Digital Surface Model (DSM)\n NEON Harvard Forest Field Site")
Data Tip: Note that when we assign break values a set of 4 values will result in 3 bins of data.
If we need to create multiple plots using the same color palette, we can create an R object (myCol) for the set of colors that we want to use. We can then quickly change the palette across all plots by simply modifying the myCol object.
We can label the x- and y-axes of our plot too using xlab and ylab.
# Assign color to a object for repeat use/ ease of changing
myCol = terrain.colors(3)
# Add axis labels
plot(DSM_HARV,
breaks = c(300, 350, 400, 450),
col = myCol,
main = "Digital Surface Model\nNEON Harvard Forest Field Site",
xlab = "UTM Westing Coordinate (m)",
ylab = "UTM Northing Coordinate (m)")
We can layer a raster on top of a hillshade raster for the same area, and use a transparency factor to created a 3-dimensional shaded effect. A hillshade is a raster that maps the shadows and texture that you would see from above when viewing terrain.
# import DSM hillshade
DSM_hill_HARV <- raster("NEON-DS-Airborne-Remote-Sensing/HARV/DSM/HARV_DSMhill.tif")
# plot hillshade using a grayscale color ramp that looks like shadows.
plot(DSM_hill_HARV,
col = grey(1:100/100), # create a color ramp of grey colors
legend = FALSE,
main = "Hillshade - DSM\n NEON Harvard Forest Field Site",
axes = FALSE)
Data Tip: Turn off, or hide, the legend on a plot
using legend=FALSE.
We can layer another raster on top of our hillshade using by using add = TRUE. Let’s overlay DSM_HARV on top of the hill_HARV.
# plot hillshade using a grayscale color ramp that looks like shadows.
plot(DSM_hill_HARV,
col = grey(1:100/100), #create a color ramp of grey colors
legend = F,
main = "DSM with Hillshade \n NEON Harvard Forest Field Site",
axes = FALSE)
# add the DSM on top of the hillshade
plot(DSM_HARV,
col = rainbow(100),
alpha = 0.4, # partly transparent
add = T, # layers this plot on top of the previous plot
legend = F)
The alpha value determines how transparent the colors will be (0 being transparent, 1 being opaque). Note that here we used the color palette rainbow() instead of terrain.color().
Original Authors: Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul
Original Last Updated: Apr 8, 2021
Sometimes we encounter raster datasets that do not “line up” when plotted or analyzed. Rasters that don’t line up are most often in different Coordinate Reference Systems (CRS).
This tutorial explains how to deal with rasters in different, known CRSs. It will walk though reprojecting rasters in R using the projectRaster() function in the raster package.
In the Plot Raster Data in R tutorial, we learned how to layer a raster file on top of a hillshade for a nice looking basemap. In that tutorial, all of our data were in the same CRS. What happens when things don’t line up?
Let’s create a map of the Harvard Forest Digital Terrain Model (DTM_HARV) draped or layered on top of the hillshade (DTM_hill_HARV).
DTM_HARV <- raster("NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_dtmCrop.tif")
DTM_hill_HARV <- raster("NEON-DS-Airborne-Remote-Sensing/HARV/DTM/HARV_DTMhill_WGS84.tif")
# plot hillshade using a grayscale color ramp
#plot(DTM_hill_HARV,
#col = grey(1:100/100),
#legend = FALSE,
#main = "DTM Hillshade\n NEON Harvard Forest Field Site")
# overlay the DTM on top of the hillshade
#plot(DTM_HARV,
#col = terrain.colors(10),
#alpha = 0.4,
#add = TRUE,
#legend = FALSE)
# I ran this chunk twice and it crashed my R session both times
# It's supposed to only plot the first one
Our results are curious - the Digital Terrain Model (DTM_HARV) did not plot on top of our hillshade. The hillshade plotted just fine on it’s own. Let’s try to plot the DTM on it’s own to make sure there are data there.
Code Tip: For boolean R elements, such as
add = TRUE, you can use T and F
in place of TRUE and FALSE
# Plot DTM
plot(DTM_HARV,
col = terrain.colors(10),
alpha = 1,
legend = F,
main = "Digital Terrain Model\n NEON Harvard Forest Field Site")
Our DTM seems to contain data and plots just fine. Let’s next check the Coordinate Reference System (CRS) and compare it to our hillshade.
# view crs for DTM
crs(DTM_HARV)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16018]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
# view crs for hillshade
crs(DTM_hill_HARV)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
Aha! DTM_HARV is in the UTM projection. DTM_hill_HARV is in Geographic WGS84 - which is represented by latitude and longitude values. Because the two rasters are in different CRSs, they don’t line up when plotted in R. We need to reproject DTM_hill_HARV into the UTM CRS. Alternatively, we could project DTM_HARV into WGS84.
We can use the projectRaster function to reproject a raster into a new CRS. Keep in mind that reprojection only works when you first have a defined CRS for the raster object that you want to reproject. It cannot be used if no CRS is defined. Lucky for us, the DTM_hill_HARV has a defined CRS.
Data Tip: When we reproject a raster, we move it from one “grid” to another. Thus, we are modifying the data! Keep this in mind as we work with raster data. To use the projectRaster function, we need to define two things:
The syntax is projectRaster(RasterObject, crs = CRSToReprojectTo)
We want the CRS of our hillshade to match the DTM_HARV raster. We can thus assign the CRS of our DTM_HARV to our hillshade within the projectRaster() function as follows: crs=crs(DTM_HARV)
# reproject to UTM
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV))
# compare attributes of DTM_hill_UTMZ18N to DTM_hill
crs(DTM_hill_UTMZ18N_HARV)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 18N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-75,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16018]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
# compare attributes of DTM_hill_UTMZ18N to DTM_hill
crs(DTM_hill_HARV)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
# compare attributes of DTM_hill_UTMZ18N to DTM_hill
extent(DTM_hill_UTMZ18N_HARV)
class : Extent
xmin : 731397.3
xmax : 733205.3
ymin : 4712403
ymax : 4713907
extent(DTM_hill_HARV)
class : Extent
xmin : -72.18192
xmax : -72.16061
ymin : 42.52941
ymax : 42.54234
Notice in the output above that the crs() of DTM_hill_UTMZ18N_HARV is now UTM. Also, the extent values of DTM_hillUTMZ18N_HARV are different from DTM_hill_HARV.
Let’s next have a look at the resolution of our reprojected hillshade.
# compare resolution
res(DTM_hill_UTMZ18N_HARV)
[1] 1.000 0.998
The output resolution of DTM_hill_UTMZ18N_HARV is 1 x 0.998. Yet, we know that the resolution for the data should be 1m x 1m. We can tell R to force our newly reprojected raster to be 1m x 1m resolution by adding a line of code (res=).
# adjust the resolution
DTM_hill_UTMZ18N_HARV <- projectRaster(DTM_hill_HARV,
crs = crs(DTM_HARV),
res = 1)
# view resolution
res(DTM_hill_UTMZ18N_HARV)
[1] 1 1
Let’s plot our newly reprojected raster.
# plot newly reprojected hillshade
plot(DTM_hill_UTMZ18N_HARV,
col = grey(1:100/100),
legend = F,
main = "DTM with Hillshade\n NEON Harvard Forest Field Site")
# overlay the DTM on top of the hillshade
plot(DTM_HARV,
col = terrain.colors(100),
alpha = 0.4,
add = T,
legend = F)
We have now successfully draped the Digital Terrain Model on top of our hillshade to produce a nice looking, textured map!
Use the
NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif
and
NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif
files.
We’ll load in the two rasterLayers from the files listed above:
DSM_SJER <- raster("NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_dsmCrop.tif")
DSM_hill_UTMZ18N_SJER <- raster("NEON-DS-Airborne-Remote-Sensing/SJER/DSM/SJER_DSMhill_WGS84.tif")
crs(DSM_SJER)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 11N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-117,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16011]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
crs(DSM_hill_UTMZ18N_SJER)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
DSM_hill_SJER_UTM <- projectRaster(DSM_hill_UTMZ18N_SJER,
crs = crs(DSM_SJER))
crs(DSM_hill_SJER_UTM)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 11N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-117,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16011]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Great work! You have learned a lot about working with raster data and plotting using base R.